#!/usr/bin/perl -w
use strict;
use DBI;

# Jason Weirather 2010.
# This script will allow you to associate gene names with hits when the gene is located within a certain distance of the hit.
# This distance is set in the upstream and downstream parameters.
# As of now you need to hard code in which genome you are using under species.

if(scalar(@ARGV) < 2) { die "cat list_of_hits.txt | ./annotate_matches_ctcf.pl genome cell_type\n"; }
my $species = shift(@ARGV);
my $celltype = shift(@ARGV);
my $upstream = 500000000;
my $downstream = 500000000;
while(my $line = <STDIN>) {
  if($line =~ /^###/) { } else {
    chomp($line);
    my @fields = split(/\t/,$line);
    find_features($fields[0],$fields[1],$fields[2],$species,$celltype,$line,$upstream,$downstream);  
  }
}

sub find_features {
  my $chromosome = shift @_;
  my $match_start = shift @_;
  my $match_finish = shift @_;
  my $match_midpoint = int(($match_start+$match_finish)/2);
  my $species = shift @_;
  my $celltype = shift @_;
  my $line = shift @_;
  my $upstream_distance = shift @_;
  my $downstream_distance = shift @_;
  my $printed_something = 0;
    #this function will try and describe where it is.  We will use the 
    #mysql refGene.sql database from UCSC for this part.  need to add multi species support later 
  my $database = "compressed_genomes";
  my $platform = "mysql";
  my $user = "x";
  my $dbh = DBI->connect("dbi:$platform:$database",$user,'');
          #is it upstream from a gene?
  ## first find the bounds.
   
  my $query = "SELECT zebraphish.ctcf_bounds.chromosome, zebraphish.ctcf_bounds.start, zebraphish.ctcf_bounds.finish FROM zebraphish.ctcf_bounds WHERE species = '$species' AND cell_type = '$celltype' AND start <= $match_midpoint AND finish >= $match_midpoint AND chromosome = '$chromosome'";
  #die "$query\n";
  my $st = $dbh->prepare($query) or die "preparing query failed: $DBI::errstr";
  $st->execute() or die "The execution of query field: $DBI::errstr";
  my $ctcf_start = -1;
  my $ctcf_finish = -1;
  while(my $row = $st->fetchrow_hashref()) {
    $ctcf_start = $row->{'start'};
    $ctcf_finish = $row->{'finish'};
  }
  #die "$ctcf_start\t$ctcf_finish\n";
  if($ctcf_start < 0 || $ctcf_finish < 0) { 
    #### we don't have any ctcf information for this chromosome so skip it.
 } else {

  $query = "SELECT strand, txStart, txEnd, name, name2, 'upstream' as relative_location FROM $species" . "_genes WHERE chrom = \"$chromosome\" AND (";
  $query .= "(strand = \"+\" AND $match_start < txStart AND $match_start > (txStart-$upstream_distance)) OR ";
  $query .= "(strand = \"-\" AND $match_finish > txEnd AND $match_finish < (txEnd+$upstream_distance))";
  $query .= ") AND (txEnd > $ctcf_start AND txStart < $ctcf_finish)";
        #echo "upstream query: $query<br>\n";
  $st = $dbh->prepare($query) or die "preparing query failed: $DBI::errstr";
  $st->execute() or die "The execution of query field: $DBI::errstr";

  $query = "SELECT strand, txStart, txEnd, name, name2, 'upstream' as relative_location FROM $species" . "_genes WHERE chrom = \"$chromosome\" AND (";
  $query .= "(strand = \"+\" AND $match_start < txStart AND $match_start > (txStart-$upstream_distance)) OR ";
  $query .= "(strand = \"-\" AND $match_finish > txEnd AND $match_finish < (txEnd+$upstream_distance))";
  $query .= ") AND (txEnd > $ctcf_start AND txStart < $ctcf_finish)";
        #echo "upstream query: $query<br>\n";
  $st = $dbh->prepare($query) or die "preparing query failed: $DBI::errstr";
  $st->execute() or die "The execution of query field: $DBI::errstr";
  my $upstream_line = "";
  my $name;
  my $name2;
  my $distance;
  while(my $row = $st->fetchrow_hashref()) {
    #see all the genes it is upstream of.
    $distance = 0;
    if($row->{'strand'} eq '+') { $distance = $row->{'txStart'}-$match_finish; }
    else { $distance = $match_start-$row->{'txEnd'}; }
    $name = $row->{'name'};
    $name2 = $row->{'name2'};
    #echo "falls $distance bp upstream of $name, $name2<br>\n";
    $upstream_line .= "$name2 ($name) / $distance bp<br>";
    print "$name2\t$name\tupstream\t$distance\t$line\n";
    $printed_something++;
  }
       ### preg_match('/^(.*)<br>$/',$upstream_line,$upmatch); # what does this do??
       ## $upstream_line = "<td>" . $upmatch[1] . "</td>"; ##
  #print "upstream\n$upstream_line\n\n";
       #is it downstream from a gene?
  $query = "SELECT strand, txStart, txEnd, name, name2, 'downstream' as relative_location FROM $species"."_genes WHERE chrom = \"$chromosome\" AND (";
  $query .= "(strand = \"-\" AND $match_start < txStart AND $match_start > (txStart-$upstream_distance)) OR ";
  $query .= "(strand = \"+\" AND $match_finish > txEnd AND $match_finish < (txEnd+$upstream_distance))";
  $query .= ") AND (txEnd > $ctcf_start AND txStart < $ctcf_finish)";
        #echo "downstream query: $query<br>\n";
  $st = $dbh->prepare($query) or die "preparing query failed: $DBI::errstr";
  $st->execute() or die "The execution of query field: $DBI::errstr";
  my $downstream_line = "";
  while(my $row = $st->fetchrow_hashref()) {
    #see all the genes it is downstream of.
    $distance = 0;
    if($row->{'strand'} eq '+') { $distance = $match_start-$row->{'txEnd'}; }
    else { $distance = $row->{'txStart'}-$match_finish; }
    $name = $row->{'name'};
    $name2 = $row->{'name2'};
          #echo "falls $distance bp downstream of $name, $name2<br>\n";
    #$downstream_line .= "$name2 ($name) / $distance bp<br>";
    print "$name2\t$name\tdownstream\t$distance\t$line\n";
    $printed_something++;
  } 
      ##  preg_match('/^(.*)<br>$/',$downstream_line,$downmatch);
      ##  $downstream_line = "<td>" . $downmatch[1] . "</td>";
      ##  fwrite($fh,$downstream_line);
  #print "downstream\n$downstream_line\n\n";
          #is it in a gene?
  $query = "SELECT strand, name, name2, 'inside' as relative_location, exonCount, exonStarts, exonEnds FROM $species"."_genes WHERE chrom = \"$chromosome\" AND (";
  $query .= "($match_finish > txStart AND $match_start < txEnd)";
  $query .= ") AND (txEnd > $ctcf_start AND txStart < $ctcf_finish)";
        #echo "inside query: $query<br>\n";
  $st = $dbh->prepare($query) or die "preparing query failed: $DBI::errstr";
  $st->execute() or die "The execution of query field: $DBI::errstr";
  my $inside_line = "";
        #while($row = mysql_fetch_object($result)) {
        #  $name = $row->name;
        #  $name2 = $row->name2;
        #  echo "falls inside $name, $name2<br>\n";
        #}
        #if it is in a gene, we need to figure out whether or not its in an exon for each gene that its in.
  while(my $row = $st->fetchrow_hashref()) {
    $name = $row->{'name'};
    $name2 = $row->{'name2'};
    my $exonCount = $row->{'exonCount'};
    my @exonStarts = split(/,/,$row->{'exonStarts'});
    my @exonEnds = split(/,/,$row->{'exonEnds'});
    $inside_line .= "$name2 ($name)";
    for(my $i = 0; $i < $exonCount; $i++) {
      #now see if it dips into any exons.
            #echo "check exons<br>\n";
      if(($match_start < $exonStarts[$i] && $match_finish > $exonStarts[$i]) || ($match_finish > $exonEnds[$i] && $match_start < $exonEnds[$i]) || ($match_start > $exonStarts[$i] && $match_finish < $exonEnds[$i]) || ($match_start < $exonStarts[$i] && $match_finish > $exonEnds[$i])) {
               #it dips into an exon
               #echo "in exon " . ($i+1) . "<br>\n";
        $inside_line .= " / exon " . ($i+1);
    print "$name2\t$name\texon\t" . ($i+1) . "\t$line\n";
    $printed_something++;
      }
            #if there is a next exon, check the intron
      if($i < $exonCount-1) {
              #see if it dips into any introns
        if(($match_start < $exonEnds[$i] && $match_finish > $exonEnds[$i]) || ($match_start < $exonStarts[$i+1] && $match_finish > $exonStarts[$i+1]) || ($match_start > $exonEnds[$i] && $match_finish < $exonStarts[$i+1]) || ($match_start < $exonEnds[$i] && $match_finish > $exonStarts[$i+1])) {
          #it dips into an intron
          #echo "in intron " . ($i+1) . "<br>\n";
          $inside_line .= " / intron " . ($i+1);
    print "$name2\t$name\tintron\t" . ($i+1) . "\t$line\n";
    $printed_something++;
        }
      }
    }
    $inside_line .= "<br>";
  }
        ##preg_match('/^(.*)<br>$/',$inside_line,$inmatch);
        ##$inside_line = "<td>" . $inmatch[1] . "</td></tr>";
        #echo "inside: $inside_line<br>\n";
        ##fwrite($fh,$inside_line);
  #print "inside\n$inside_line\n\n";
  if($printed_something == 0) { print "\t\t\t\t$line\n"; }
 }
  return;
}
